#include "mt3d.h"

void mt3d::solve()
{
    // Preconditioning the iteration method
    /*
    std::clog << "Constructing the preconditioner ... ";

    clcg_incomplete_Cholesky(A_PDE_, L_); 
    LT_ = L_.transpose();
    //incomplete_LU(A_PDE_, L_, LT_);

    std::clog << "done\n";
    */
    /*
    for (size_t i = 0; i < edge_num_; i++)
    {
        for (Eigen::SparseMatrix<complex_d, Eigen::RowMajor>::InnerIterator it(L_, i); it; ++it) // 列循环
        {
            if (it.value() != _zero)
            {
                std::cout << it.row() << " " << it.col() << " " << it.value() << std::endl;
            }
        }
    }
    */

    
    // Initialize the solution vector using \diag(A)
    
    P_.resize(edge_num_);
    for (size_t i = 0; i < edge_num_; i++)
    {
        for (Eigen::SparseMatrix<complex_d, Eigen::RowMajor>::InnerIterator it(A_PDE_, i); it; ++it)
        {
            if (it.row() == it.col()) P_[i] = _one/it.value();
        }
    }
    /*
    X_TE_.resize(edge_num_); X_TE_.setZero();
    X_TM_.resize(edge_num_); X_TM_.setZero();

    for (size_t i = 0; i < edge_num_; i++)
    {
        X_TE_[i] = B_TE_[i] * P_[i];
        X_TM_[i] = B_TM_[i] * P_[i];
    }
    */

    //P_.resize(edge_num_);
    X_TE_.resize(edge_num_); X_TE_.setZero();
    X_TM_.resize(edge_num_); X_TM_.setZero();
    // Call conjugate gradient to solve the linear system
    clcg_para my_para = clcg_default_parameters();
    //my_para.abs_diff = 1;
    my_para.epsilon = 1e-8;
    set_clcg_parameter(my_para);
    set_report_interval(100);

    std::clog << "Solving the PDE in TE mode ... \n";
    MinimizePreconditioned(X_TE_, B_TE_, CLCG_PCG);

    std::clog << "Solving the PDE in TM mode ... \n";
    MinimizePreconditioned(X_TM_, B_TM_, CLCG_PCG);

    return;
}
